home *** CD-ROM | disk | FTP | other *** search
/ AP Professional Graphics CD-ROM Library / AP Professional Graphics CD-ROM Library.iso / pc / unix / appendix / gemsiii / partit3d / prtition.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-03-17  |  12.1 KB  |  373 lines
  1. /* partition.c: module to partition 3D convex face with an arbitrary plane.
  2.  * Copyright (c) Norman Chin
  3.  */ 
  4. #include <stdio.h>        /* NULL */
  5. #include <assert.h>        /* assert() */
  6. #include <malloc.h>        /* malloc() */
  7. #include <math.h>        /* fabs() */
  8. #include "GraphicsGems.h"    
  9. #include "partition.h"
  10.  
  11. /* local functions */
  12. static VERTEX *findNextIntersection(/* VERTEX *vstart, double aa,
  13.                        double bb, double cc, double dd,
  14.                        double *ixx, double *iyy, double *izz,
  15.                        SIGN *sign */);
  16. static SIGN anyEdgeIntersectWithPlane(/* double x1, double y1, double z1,
  17.                      double x2, double y2, double z2,
  18.                      double aa, double bb, double cc,
  19.                      double dd, 
  20.                      double *ixx, double *iyy, double *izz
  21.                        */);
  22. static FACE *createOtherFace(/* FACE *face, VERTEX *v1, double ixx1, 
  23.                 double iyy1, double izz1, VERTEX *v2,
  24.                 double ixx2, double iyy2, double izz2 */);
  25. static SIGN whichSideIsFaceWRTplane(/* FACE *face, double aa, double bb,
  26.                        double cc, double dd */);
  27. static VERTEX *allocVertex(/* double xx, double yy, double zz */);
  28. static FACE *allocFace(/* FACE *face, VERTEX *vhead */);
  29.  
  30. /* Partitions a 3D convex polygon (face) with an arbitrary plane into its 
  31.  * negative and positive fragments, if any, w.r.t. the partitioning plane.
  32.  * Note that inputFace is unusable afterwards since its vertex list has been
  33.  * parceled out to the other faces. It's set to null to avoid dangling
  34.  * pointer problem. 
  35.  */
  36. void partitionFaceWithPlane(aa,bb,cc,dd,inputFace,faceOn,faceNeg,facePos)
  37. double aa,bb,cc,dd;        /* partitioning plane's coefficients */
  38. FACE **inputFace;        /* face to be partitioned */
  39. FACE **faceOn;            /* returns face embedded in plane, if any  */
  40. FACE **faceNeg;            /* returns face on negative side, if any */
  41. FACE **facePos;            /* returns face on positive side, if any */
  42. {
  43.    VERTEX *v1, *v2;
  44.    double ixx1,iyy1,izz1, ixx2,iyy2,izz2;
  45.    SIGN signV1, signV2;
  46.  
  47.    *faceOn= *faceNeg= *facePos= NULL_FACE;
  48.  
  49.    /* find first intersection */
  50.    v1= findNextIntersection((*inputFace)->vhead,aa,bb,cc,dd, 
  51.                 &ixx1,&iyy1,&izz1,&signV1);
  52.    if (v1 != NULL_VERTEX) {
  53.       assert(signV1 != ZERO);
  54.  
  55.       /* first one found, find the 2nd one, if any */
  56.       v2= findNextIntersection(v1->vnext,aa,bb,cc,dd, 
  57.                    &ixx2,&iyy2,&izz2,&signV2);
  58.  
  59.       /* Due to numerical instabilities, both intersection points may
  60.        * have the same sign such as in the case when splitting very close
  61.        * to a vertex. This should not count as a split.
  62.        */
  63.       if (v2 != NULL_VERTEX && signV1 == signV2) v2= NULL_VERTEX; 
  64.  
  65.    }
  66.    else v2= NULL_VERTEX;    /* No first intersection found,
  67.                                  * therefore no second intersection.
  68.                  */
  69.                                 
  70.    /* an intersection? */
  71.    if (v1 != NULL_VERTEX && v2 != NULL_VERTEX) { /* yes, intersection */
  72.       FACE *newOtherFace;
  73.  
  74.       /* inputFace's vertex list will be modified */
  75.       newOtherFace= createOtherFace(*inputFace,v1,ixx1,iyy1,izz1,
  76.                                v2,ixx2,iyy2,izz2);
  77.       /* return split faces on appropriate lists */
  78.       if (signV1 == NEGATIVE) { 
  79.      *faceNeg= *inputFace; 
  80.          *facePos= newOtherFace;
  81.       }
  82.       else {
  83.      assert(signV1 == POSITIVE);
  84.      *faceNeg= newOtherFace;
  85.      *facePos= *inputFace;
  86.       }
  87.                     
  88.    }
  89.    else {            /* no intersection  */
  90.       SIGN side;
  91.  
  92.       /* Face is embedded or wholly to one side of partitioning plane. */
  93.       side= whichSideIsFaceWRTplane(*inputFace,aa,bb,cc,dd);
  94.       if (side == NEGATIVE)
  95.      *faceNeg= *inputFace;
  96.       else if (side == POSITIVE)
  97.      *facePos= *inputFace;
  98.       else {
  99.      assert(side == ZERO);
  100.      *faceOn= *inputFace;
  101.       }
  102.    }
  103.    /* inputFace's vertex list has been parceled out to other lists so
  104.     * set this to null.
  105.     */
  106.    *inputFace= NULL_FACE;        
  107.  
  108.    /* If a face is embedded in plane then there must not be any faces on
  109.     *    either side.
  110.     * Otherwise the face isn't embedded in plane so there must be at least 
  111.     *    one face on either side or both.
  112.     */
  113.    assert( (*faceOn != NULL_FACE) ? 
  114.         (*faceNeg == NULL_FACE &&
  115.          *facePos == NULL_FACE) :
  116.         (*faceNeg != NULL_FACE ||
  117.          *facePos != NULL_FACE) );
  118. } /* partitionFaceWithPlane() */
  119.                 
  120. /* Finds next intersection on or after vstart. 
  121.  * 
  122.  * If an intersection is found, 
  123.  *    a pointer to first vertex of the edge is returned, 
  124.  *    the intersection point (ixx,iyy,izz) and its sign is updated. 
  125.  * Otherwise a null pointer is returned.
  126.  */
  127. static VERTEX *findNextIntersection(vstart,aa,bb,cc,dd,ixx,iyy,izz,sign)
  128. VERTEX *vstart;            /* where to start searching on vertex list */
  129. double aa,bb,cc,dd;        /* plane's coefficients */
  130. double *ixx,*iyy,*izz;        /* returned intersection point */
  131. SIGN *sign;        /* returned classification of intersection point */
  132. {
  133.    VERTEX *vtrav;
  134.  
  135.    /* for all edges starting from vstart ... */
  136.    for (vtrav= vstart; vtrav->vnext != NULL_VERTEX; vtrav= vtrav->vnext) {
  137.       if ((*sign= anyEdgeIntersectWithPlane(vtrav->xx,vtrav->yy,vtrav->zz,
  138.                         vtrav->vnext->xx,vtrav->vnext->yy,
  139.                         vtrav->vnext->zz,aa,bb,cc,dd,
  140.                         ixx,iyy,izz))) {
  141.      return(vtrav);
  142.       }
  143.    }
  144.  
  145.    return(NULL_VERTEX);
  146. } /* findNextIntersection() */
  147.  
  148. /* Memory allocated for split face's vertices and pointers tediously updated.
  149.  */
  150. static FACE *createOtherFace(face,v1,ixx1,iyy1,izz1,v2,ixx2,iyy2,izz2)
  151. FACE *face;            /* face to be split */
  152. VERTEX *v1;    /* 1st vertex of edge of where 1st intersection was found */
  153. double ixx1,iyy1,izz1;        /* 1st intersection  */
  154. VERTEX *v2;     /* 1st vertex of edge of where 2nd intersection was found */
  155. double ixx2,iyy2,izz2;        /* 2nd intersection */
  156. {
  157.    VERTEX *i1p1, *i2p1;        /* new vertices for original face  */
  158.    VERTEX *i1p2, *i2p2;        /* new vertices for new face */
  159.    VERTEX *p2end;        /* new vertex for end of new face */
  160.    VERTEX *vtemp;        /* place holders */
  161.    register VERTEX *beforeV2;    /* place holders */
  162.    FACE *newFace;        /* newly allocated face */
  163.  
  164.    /* new intersection vertices */
  165.    i1p1= allocVertex(ixx1,iyy1,izz1); 
  166.    i2p1= allocVertex(ixx2,iyy2,izz2);
  167.    i1p2= allocVertex(ixx1,iyy1,izz1);
  168.    i2p2= allocVertex(ixx2,iyy2,izz2); 
  169.  
  170.    /* duplicate 1st vertex of 2nd list to close it up */
  171.    p2end= allocVertex(v2->xx,v2->yy,v2->zz);
  172.  
  173.    vtemp= v1->vnext;
  174.  
  175.    /* merge both intersection vertices i1p1 & i2p1 into 1st list */
  176.    if (ISVERTEX_EQ(i2p1,v2->vnext)) { /* intersection vertex coincident? */
  177.       FREEVERTEX(i2p1);        /* yes, so free it */
  178.       i1p1->vnext= v2->vnext;
  179.    }
  180.    else {
  181.       i2p1->vnext= v2->vnext;    /* attach intersection list onto 1st list */
  182.       i1p1->vnext= i2p1;    /* attach both intersection vertices */
  183.    }
  184.    v1->vnext= i1p1; /* attach front of 1st list to intersection vertices */
  185.  
  186.    /* merge intersection vertices i1p2, i2p2 & p2end into second list */
  187.    i2p2->vnext= i1p2;        /* attach both intersection vertices */
  188.    v2->vnext= i2p2;        /* attach 2nd list to interection vertices */
  189.    if (vtemp == v2) {
  190.       i1p2->vnext= p2end;    /* close up 2nd list */
  191.    }
  192.    else {
  193.       if (ISVERTEX_EQ(i1p2,vtemp)) { /* intersection vertex coincident? */
  194.      FREEVERTEX(i1p2);    /* yes, so free it */
  195.      i2p2->vnext= vtemp;    /* attach intersection vertex to 2nd list */
  196.       }
  197.       else {
  198.      i1p2->vnext= vtemp;    /* attach intersection list to 2nd list */
  199.       }
  200.       /* find previous vertex to v2 */
  201.       for (beforeV2= vtemp; beforeV2->vnext != v2; beforeV2= beforeV2->vnext)
  202.      ;            /* lone semi-colon */
  203.       beforeV2->vnext= p2end;    /* and attach it to complete the 2nd list */
  204.    }
  205.  
  206.    /* copy original face info but with new vertex list */
  207.    newFace= allocFace(face,v2);
  208.  
  209.    return(newFace);
  210. } /* createOtherFace() */
  211.  
  212. /* Determines which side a face is with respect to a plane with coefficient
  213.  * (aa,bb,cc,dd).
  214.  *
  215.  * However, due to numerical problems, when a face is very close to the plane,
  216.  * some vertices may be misclassified. 
  217.  * There are several solutions, two of which are mentioned here:
  218.  *    1- classify the one vertex furthest away from the plane, (note that
  219.  *       one need not compute the actual distance) and use that side.
  220.  *    2- count how many vertices lie on either side and pick the side
  221.  *       with the maximum. (this is the one implemented).
  222.  */
  223. static SIGN whichSideIsFaceWRTplane(face,aa,bb,cc,dd)
  224. FACE *face;
  225. double aa,bb,cc,dd;
  226. {
  227.    register VERTEX *vtrav;
  228.    double value;
  229.    boolean isNeg, isPos;
  230.  
  231.    isNeg= isPos= FALSE;
  232.    
  233.    for (vtrav= face->vhead; vtrav->vnext != NULL_VERTEX; vtrav= vtrav->vnext){
  234.       value= (aa*vtrav->xx) + (bb*vtrav->yy) + (cc*vtrav->zz) + dd;
  235.       if (value < -TOLER) 
  236.      isNeg= TRUE;
  237.       else if (value > TOLER)
  238.      isPos= TRUE;
  239.       else assert(ISDOUBLE_EQ(value,0.0));
  240.    } /* for vtrav */ 
  241.  
  242.    /* in the very rare case that some vertices slipped thru to other side of
  243.     * plane due to round-off errors, execute the above again but count the 
  244.     * vertices on each side instead and pick the maximum.
  245.     */
  246.    if (isNeg && isPos) {    /* yes so handle this numerical problem */
  247.       int countNeg, countPos;
  248.       
  249.       /* count how many vertices are on either side */
  250.       countNeg= countPos= 0;
  251.       for (vtrav= face->vhead; vtrav->vnext != NULL_VERTEX; 
  252.        vtrav= vtrav->vnext) {
  253.      value= (aa*vtrav->xx) + (bb*vtrav->yy) + (cc*vtrav->zz) + dd;
  254.      if (value < -TOLER)
  255.         countNeg++;
  256.      else if (value > TOLER)
  257.         countPos++;
  258.      else assert(ISDOUBLE_EQ(value,0.0));
  259.       } /* for */
  260.  
  261.       /* return the side corresponding to the maximum */
  262.       if (countNeg > countPos) return(NEGATIVE);
  263.       else if (countPos > countNeg) return(POSITIVE);
  264.       else return(ZERO);
  265.    }
  266.    else {            /* this is the usual case */
  267.       if (isNeg) return(NEGATIVE);
  268.       else if (isPos) return(POSITIVE);
  269.       else return(ZERO);
  270.    }
  271. } /* whichSideIsFaceWRTplane() */
  272.  
  273. /* Determines if an edge bounded by (x1,y1,z1)->(x2,y2,z2) intersects
  274.  * the plane with the coefficients (aa,bb,cc,dd).
  275.  * 
  276.  * If there's an intersection, 
  277.  *    the sign of (x1,y1,z1), NEGATIVE or POSITIVE, w.r.t. the plane is
  278.  *    returned with the intersection (ixx,iyy,izz) updated.
  279.  * Otherwise ZERO is returned.    
  280.  */
  281. static SIGN anyEdgeIntersectWithPlane(x1,y1,z1,x2,y2,z2,aa,bb,cc,dd,
  282.                     ixx,iyy,izz)
  283. double x1,y1,z1, x2,y2,z2;    /* both vertices of edge */
  284. double aa,bb,cc,dd;        /* partitioning plane's coefficients */
  285. double *ixx,*iyy,*izz;        /* returned intersection point, if any */
  286. {
  287.    double temp1, temp2;
  288.    int sign1, sign2;        /* must be int since gonna do a bitwise ^ */
  289.  
  290.    /* get signs */
  291.    temp1= (aa*x1) + (bb*y1) + (cc*z1) + dd;
  292.    if (temp1 < -TOLER)
  293.       sign1= -1;
  294.    else if (temp1 > TOLER)
  295.       sign1= 1;
  296.    else {
  297.       /* edges beginning with a 0 sign are not considered valid intersections
  298.        * case 1 & 6 & 7, see text
  299.        */
  300.       assert(ISDOUBLE_EQ(temp1,0.0));
  301.       return(ZERO);
  302.    }
  303.  
  304.    temp2= (aa*x2) + (bb*y2) + (cc*z2) + dd;
  305.    if (temp2 < -TOLER)
  306.       sign2= -1;
  307.    else if (temp2 > TOLER)
  308.       sign2= 1;
  309.    else {            /* case 8 & 9, see text  */
  310.       assert(ISDOUBLE_EQ(temp2,0.0));
  311.       *ixx= x2;
  312.       *iyy= y2;
  313.       *izz= z2;
  314.  
  315.       return( (sign1 == -1) ? NEGATIVE : POSITIVE);
  316.    }
  317.  
  318.    /* signs different? 
  319.     * recall: -1^1 == 1^-1 ==> 1    case 4 & 5, see text
  320.     *         -1^-1 == 1^1 ==> 0    case 2 & 3, see text
  321.     */
  322.    if (sign1 ^ sign2) {
  323.       double dx,dy,dz;
  324.       double denom, tt;
  325.  
  326.       /* compute intersection point */
  327.       dx= x2-x1;
  328.       dy= y2-y1;
  329.       dz= z2-z1;
  330.  
  331.       denom= (aa*dx) + (bb*dy) + (cc*dz);
  332.       tt= - ((aa*x1) + (bb*y1) + (cc*z1) + dd) / denom;
  333.  
  334.       *ixx= x1 + (tt * dx);
  335.       *iyy= y1 + (tt * dy);
  336.       *izz= z1 + (tt * dz);
  337.  
  338.       assert(sign1 == 1 || sign1 == -1);
  339.  
  340.       return( (sign1 == -1) ? NEGATIVE : POSITIVE );
  341.    }
  342.    else return(ZERO);
  343.  
  344. } /* anyEdgeIntersectWithPlane() */
  345.  
  346. /* allocates vertex with coordinates (xx,yy,zz) */
  347. static VERTEX *allocVertex(xx,yy,zz)
  348. double xx,yy,zz;
  349. {
  350.    VERTEX *newVertex;
  351.  
  352.    newVertex= NEWTYPE(VERTEX); assert(newVertex != NULL_VERTEX);
  353.    newVertex->xx= xx; newVertex->yy= yy; newVertex->zz= zz; 
  354.    newVertex->vnext= NULL_VERTEX; 
  355.    
  356.    return(newVertex);
  357. } /* allocVertex() */
  358.  
  359. /* allocates face with list of vertices, vhead */
  360. static FACE *allocFace(face,vhead)
  361. FACE *face;
  362. VERTEX *vhead;
  363. {
  364.    FACE *newFace;
  365.  
  366.    newFace= NEWTYPE(FACE); assert(newFace != NULL_FACE); 
  367.    *newFace= *face;
  368.    newFace->vhead= vhead;
  369.  
  370.    return(newFace);
  371. } /* allocFace() */
  372.  
  373.